# Dr. M. Baron, Statistical Machine Learning class, STAT-427/627

# DIMENSION REDUCTION: 
# PRINCIPAL COMPONENTS AND PARTIAL LEAST SQUARES

# 1. Principal Components Regression

# Import necessary libraries
! pip install pandas;
! pip install numpy;
! pip install scikit-learn;
! pip install matplotlib;
! pip install ISLP;

from ISLP import load_data
import numpy as np
import pandas as pd
from sklearn.model_selection import train_test_split, cross_val_predict
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.cross_decomposition import PLSRegression
from sklearn.metrics import mean_squared_error
import matplotlib.pyplot as plt
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score

# Prepare the Auto dataset 
auto = load_data('auto')

# Remove the 'name' column and encode 'origin' as categorical
auto = auto.drop(columns=['name'])
auto['origin'] = auto['origin'].astype('category')

# Create dummy variables for 'origin'
auto = pd.get_dummies(auto, drop_first=True)

# Define X (predictors) and y (response)
X = auto.drop(columns=['mpg'])
y = auto['mpg']

# Standardize the predictors
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)

# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y, test_size=0.5, random_state=42)

### Principal Components Regression (PCR) ###

# Fit PCA on the training data
pca = PCA()
X_train_pca = pca.fit_transform(X_train)

# Fit PCR (regression on principal components)
reg = LinearRegression()
reg.fit(X_train_pca, y_train)

# Number of components
n_components = X_train_pca.shape[1]

# Variance explained by each component
explained_variance = pca.explained_variance_ratio_

# Cumulative variance explained
cumulative_variance = np.cumsum(explained_variance)

# Print the variance explained
print("Explained variance by each component:", explained_variance)
print("Cumulative variance:", cumulative_variance)

# Scree plot (variance explained by principal components)
plt.figure(figsize=(8, 6))
plt.plot(np.arange(1, n_components + 1), cumulative_variance, marker='o')
plt.title('Scree Plot')
plt.xlabel('Number of Principal Components')
plt.ylabel('Cumulative Variance Explained')
plt.show()

# Cross-validation for PCR
y_pred_cv = cross_val_predict(LinearRegression(), X_train_pca, y_train, cv=10)
mse_cv = mean_squared_error(y_train, y_pred_cv)
print("PCR Cross-validated MSE:", mse_cv)

Explained variance by each component: [0.57687438 0.15195646 0.11176798 0.08183767 0.04870938 0.01577252
 0.00887124 0.00421037]
Cumulative variance: [0.57687438 0.72883084 0.84059882 0.92243649 0.97114587 0.98691839
 0.99578963 1.        ]

PCR Cross-validated MSE: 12.954212115796144

### Partial Least Squares (PLS) ###

# Fit PLS with cross-validation to determine the optimal number of components
pls = PLSRegression(n_components=6)
pls.fit(X_train, y_train)

# Predict on test data
y_pred_pls = pls.predict(X_test)
pls_mse = mean_squared_error(y_test, y_pred_pls)
print("PLS Test MSE:", pls_mse)

# Cross-validated PLS
y_pred_pls_cv = cross_val_predict(pls, X_train, y_train, cv=10)
pls_cv_mse = mean_squared_error(y_train, y_pred_pls_cv)
print("PLS Cross-validated MSE:", pls_cv_mse)

# Plot cross-validated errors for PLS
plt.figure(figsize=(8, 6))
plt.plot(np.arange(1, len(y_pred_pls_cv) + 1), np.sort(np.abs(y_train - y_pred_pls_cv)))
plt.title('PLS Cross-validated Errors')
plt.xlabel('Data Points')
plt.ylabel('Absolute Prediction Error')
plt.show()

PLS Test MSE: 11.639369181139198
PLS Cross-validated MSE: 12.970407248015107